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This paper presents a probabilistic one- dimensional finite element model 
for heat transfer processes in porous heat exchangers. The Galerkin approach is 
used to develop the finite element matrices. Some of the submatrices are asymmet- 
ric due to the presence of the flow term. The Neumann expansion is used to write 
the temperature distribution as a series of random variables, and the expectation 
operator is applied to obtain the mean and deviation statistics. To demonstrate the 
feasibility of the formulation, a one-dimensional model of heat transfer phenomenon 
in superfluid flow through a porous media is considered. Results of this formulation 
agree well with the Monte-Carlo simulations and the analytical solutions. Although 
the numerical experiments are confined to parametric random variables, a formula- 
tion is presented to account for the random spatial variations. 


INTRODUCTION 

Porous heat exchangers are key components in many engineering systems 
such as high performance regenerative heat exchangers, thermal energy storage sys- 
tems, cryocoolers, and packed beds. Several techniques have been developed to 
analyze these systems. These techniques include, among others, analytical tech- 
niques, such as separation of variables [1,2], Riemann method [3] and similarity 
transformation [4]; -semi-analytical techniques such as orthogonal collocation [5] 
and collocation-perturbation [6,7]; and numerical methods such as numerical inte- 
gration [8], shooting and Runge-Kutta integration [9], and finite element methods 
[10]. In these techniques the above systems are considered as deterministic. That 
is, the problems are formulated in terms of mean- values of the properties neglecting 
variations in the mean values. Experimental measurements, however, show that 
the properties of the systems may vary significantly in a random fashion, especially 
near a low temperature. Given the stringent demand on the design of modern heat 
exchangers, these deterministic models may not be adequate. 

Random properties can be incorporated in the above techniques using the prob- 
ability theories and the theories of differential equations [11-14]. In many applica- 
tions, it is difficult to solve the resulting differential equations in closed form even 
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when the randomness is not considered [15]. In the last 30 years. Finite Element 
Methods (FEMs) have successfully been applied to solve a large class of determinis- 
tic systems, and more recently, probabilistic systems. Current trends for analyzing 
random systems in engineering are given by Iyengar and Dash [16], Ibrahim [17], 
Shinozuka [18], and Benaroya and Rehak [19]. Chamis and coworkers [20, 21] have 
developed a general purpose finite element computer code called PICAN for prob- 
abilistic analysis of composite materials. 

The numerical techniques available to solve problems consisting of random vari- 
ables and functions may be broadly classified into two categories; statistical and 
nonstatistical. Most statistical techniques rely on numerical simulations among 
which Monte Carlo simulation has been widely used [22]. The nonstatistical tech- 
niques include perturbation methods [23-26], spectral decomposition methods [15, 
27], and basis random variable methods [28, 29]. 

In addition to the FEMs, several investigators have used the Boundary Element 
Methods (BEMs). The random operator problems were analyzed by Ettouney et 
al. [30-32], and more recently by Manolis and Shaw [33]. Burczynski [34] employed 
the direct boundary element method to develop two distinct procedures for the 
treatment of random potential problems. Cheng and Lafe [35] employed the indi- 
rect boundary element method to obtain stochastic integral equations for boundary 
potentials and fluxes in terms of fictitious boundary sources. Other applications of 
the boundary element method include the first order perturbation method devel- 
oped by Drewniak [36] for the analysis of heat conduction problems with random 
heat transfer and random heat conduction coefficient, respectively, and a procedure 
for the analysis of time dependent problems in the frequency domain developed by 
Burczynski and John [37]. 

The probabilistic methods discussed here, however, are largely confined to struc- 
tural systems, and very little effort has been made to develop methods for porous 
heat exchangers. In this paper, a probabilistic one-dimensional finite element model 
for heat transfer process in porous heat exchangers is presented. The formulation 
is based on the Galerkin method, the spectral decomposition of random processes, 
and the Neumann expansion. 


MATHEMATICAL FORMULATION 

This section is divided into three parts: deterministic finite element model, prob- 
abilistic finite element models for parametric randomness and stochastic processes, 
and Monte-Carlo models. These models are considered next. 

Deterministic Finite Element Model .- In order to develop a deterministic 
analytical model for the heat transfer process in fluid flow through porous media, 
consider the schematic of a one dimensional heat exchanger as shown in Figure 1. 
It is assumed that the system is stationary; that is. the system parameters are not 
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changing with time. Using energy balance, one can derive the following differential 
equations 


. cPT . tjL dT 

kf — + hs(t - T) = c— 


(1) 


k s — + fis{T-t)=0 (2) 

where kj and k s are the thermal conductivity coefficients of the fluid and the solid, 
h is the convective heat transfer coefficient, s is the porosity coefficient, T and t 
are the temperature distributions in the fluid and the solid, and the coefficient c 
accounts for the energy transfer due to fluid flow. In addition to these equations, the 
boundary conditions are also required. The formulation presented here is applicable 
to various boundary conditions. However, for simplicity, it is assumed that the 
terminal temperatures of the fluid and the solid at the two ends are prescribed. 
These conditions are written as 


T(x = 0) = To, t(x = 0) = t o T(z = T) = T l , and t(x = L) =t L (3) 

where L is the length of the heat exchanger which is discretized into several finite 
elements. If N is the vector of the shape functions defined over the elements, then 
the temperatures T and t may be approximated as 

T = N T T n and t = N T t n (4) 


where T„ and t n represent the vector of the nodal temperatures of the fluid and 
the solid, and the superscript T represents the transpose. Observe that both solid 
and fluid regions have been discretized into an equal number of elements, and the 
same shape functions have been used for both temperature distributions. This is 
not necessary and a formulation that considers different numbern of elements and 
different shape functions is possible. Using the Galerkin approach, equations (1) 
and (2) may be written as 



<PT 



dx = 0 


(5) 



^ + hs ( T - ( ) 


dx = 0 


( 6 ) 


Performing integration by parts on some terms of equations (5) and (6), and rear- 
ranging the terms, one obtains 


E 


Tn 


t n 



(7) 
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where E is the global coefficient matrix defined as 


E = 


—kfA — hsB — cC hsB 

hsB —k s A — hsB 


( 8 ) 


and vector b appears due to partial integration of some terms in Eq. (3) and the 
boundary conditions. Matrices A, B, and C in equation (8) are defined as 


r L ,dN ..dN . T , 

a = L ( d7 ,( d7 ) dx 

(9) 

B = [ L NN T dx 

(10) 

Jo 

C = [ N(^-) T dx 
Jo ax 

(11) 


Equation (7) provides the desired deterministic finite element model. Observe that 
matrices A and B are symmetric positive definite finite element matrices, whereas 
matrix C is an asymmetric matrix. This makes matrix E asymmetric. Therefore, 
one should not use a symmetric simultaneous equation solver to solve equation (7). 

Probabilistic Finite Element Model .- As stated earlier, two types of ran- 
dom behavior may appear in the system; parametric and spatial. These two random 
processes will be considered separately. 

Parametric Randomness .- For simplicity, only k s is considered as a random 
parameter. If other parameters also vary randomly, then the formulation can be 
extended appropriately. The random parameter k s may be written as k s = k, 0 + e. 
where k s0 is the mean value of k„ and e represents the random variations with mean 
zero and standard deviation a (i.e. < e >= 0, and < e 2 >= a 2 ). Substituting the 
expression for k s into Eq. (7), one obtains 


[E\ + eE 2 ] 


T n 

t n 


= b 


( 12 ) 


where matrix E\ is the same as matrix E in equation (8) except that k s is replaced 
by fc s o and matrix E 2 is given as 


E 2 


0 0 
0 -A 


(13) 


Observe that T n and t n are now vectors of random variables. Using the Neumann 
expansion, the temperature vector can be written as 


T n 

t n 


I - e(E^E 2 ) + € 2 {E~ x E 2 ) 2 - e 3 (£r'E 2 ) 3 + • • •] E' x b (14) 
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provided that ||e(£’ 1 ^ 2 )! I < 1 is satisfied, which is reasonable for most practical 
systems. From equation (14), the expected value of the temperature vector is 

>= [/-<€> (E^E 2 )+ <e 2 > (E^E 2 ) 2 - < e 3 > (E^E 2 ) 3 + • • •] E?b 

(15) 

where <> is the expectation operator. Observe that E l 1 E 2 is constant. Further- 
more, given the probability distribution, < e* > (i = 1,2, •••,) can be computed 
numerically and in some cases analytically. Substituting these values in equation 
(15), the expected values of the temperatures can be obtained. Similarly, the second 
order characteristics of T n and t n can be obtained as follows: 

>=< b T E~ T [l + eE~ l E 2 }~ T [l + eE~ l E 2 ]~ l E~ l b > (16) 

Once again, the matrices containing e in equation (16) can be expanded in Neumann 
series to obtain the covariance matrix for temperature distribution. 

Equations (15) and (16) provide a probabilistic model for parametric randomness 
in k a . A similar approach can be used for other random parameters. 

Spatial Randomness .- Consider that k s varies randomly from point to 
point along the length of the heat exchanger and that other properties are constant. 
k s can be written as k, = k, m + k sr , where k 3m represents the mean function and k sr 
represents the stationary Gaussian process with zero mean functions and specified 
correlation function R(x,u), which is symmetric and positive definite. Using the 
Karhunen-Loeve (KL) expansion, k ST can be represented as [38] 

k S r = y ' 

X 

where {<f>i\i = 1,2,...} is a set of orthogonal eigenfunctions of certain differential 
equations and e, (i = 1,2,...) are uncorrelated random variables. These eigenfunc- 
tions satisfy the following integral equation 

fL 

A i<f>i(x) = / R(x,u)<f>i(u)du , (18) 

Jo 

where A, is an eigenvalue associated with <j>i(x). Furthermore, the coefficients e, 
(i = 1. • • • , 00 ) satisfy the following identity 

A iSij =< e.Cj >, (19) 

where < e 2 >= A, gives measure of randomness along the <p,(x) coordinate. One of 
the advantages of the series expansion is that it provides a second moment charac- 
terization of k ST in terms of uncorrelated random variables. 


(17) 
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Using Eqs. (5), (6), and (17), and following the approach presented above, one 
obtains 

—kfA — hsB — cC hsB T n i 

' = b (20) 

hsB — A — hsB t n *■ 

where matrix D x is defined as 

_ f L , . ..dN. T/ dN. , . . 

D ' = Jo ( dT )c?X (21) 

Equation (20) is very similar to Eq. (12). Therefore, the statistical characteristics 

of the temperature vector Tj can be obtained using the Neumann expan- 

sion and the procedure discussed for parametric randomness. 

The above discussion provides a probabilistic model when k s represents a stochas- 
tic process. A similar approach can be used for other stochastic processes. 

Monte-Carlo Method .- Monte-Carlo simulations rely on equations (12) 
and (20). In this technique, a random number generator is used to obtain a large 
set of random numbers that represent the desired probability distribution curve of 
the random variables. This process is repeated for each random variable. Depend- 
ing on parametric or spatial randomness, equation (12) or (20) is used to obtain 
an equal number of sets of nodal temperatures, which are then used to obtain the 
statistics for the nodal temperatures. For an accurate answer, this scheme requires 
a large number of numerical tests. This number can be reduced using the following 
approach: (1) grouping the random data, (2) performing only one test for each 
group, and (3) using the probability information to account for other data in the 
group. This approach can significantly reduce the number of numerical runs for 
accurate results. 

NUMERICAL RESULTS AND DISCUSSIONS 

To validate the formulation developed here, a dilution refrigerator heat ex- 
changer consisting of superfiuid Helium II as fluid and the sintered copper as 
the solid was considered. The system response was obtained using this scheme 
and an analytical scheme. For numerical simulations, the following parameters 
were considered: Kj = 7X10 4 W/(m.K), K s = 500W/(m.A'), h = 1200IU/(m 2 A'), 
c = ll.7oW/(m.K), and s = 0.2792m which axe typical of this system. The value 
of L was taken as 10 m. For convenience, the non-dimensional temperatures of 
the fluid Tj and the solid t, were defined as follows: 7} = (T — Tq)/{Tl — T 0 ) and 
t 3 = (t — T 0 )/(T l — T 0 ). The following boundary conditions were taken: T/(0) = 0.0, 
t a (0) = 0.8, and Tj(L ) = t a (L) = 1.0. The values of k s and h were varied to 
study the effects of these parameters on the temperature distribution. To study the 
probabilistic effects, the system statistical response was obtained using Monte-Caxlo 
simulations, the analytical schemes, and the proposed scheme. Results of this study 
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are shown below. 

Figure 2 compares the temperature distributions obtained using the analytical 
and the proposed schemes. The two results agree very well. The fluid temperature 
changes almost linearly due to the large conductivity of superfluid Hell. The solid 
media temperature decreases faster near the inlet point than inside the exchanger 
because of the strong convection near the inlet. Since the outlet point tempera- 
ture is higher than the inlet temperature, the temperature of both fluid and solid 
increases after some distance. Once the temperature of the fluid and the solid be- 
come equal, convection stops and the temperature of both media increases at the 
same rate. 

Figure 3 shows the temperature profiles of the solid and the fluid for k, = 250, 
500, 1000, 1500, and 2000 W/(m.K). It is clear that an increase in solid thermal 
conductivity causes the temperature of the solid to increase. This is because con- 
vection becomes less dominant at higher values of k s . Due to large fluid thermal 
conductivity, the fluid temperature profile remains unchanged. 

The temperature profiles for h = 600, 1200, 2400, and 12,000 W(m 2 .K) are 
shown in figure 4. As expected, increase in convective heat transfer coefficient 
causes the solid temperature to decrease rapidly and merge with the fluid temper- 
ature sooner. 

In this investigation, parameter h was considered as a random variable, and all 
other parameters were kept the same. Using a random number generator, 20,000 
random sample points having 10 % variations of mean convective heat transfer co- 
efficient with 90 % confidence were generated. These data were used in the Monte- 
Carlo method, the exact solution, and this scheme to predict the mean response 
of the temperature profile. All three schemes gave the same results. To compare 
the relative accuracy of the current scheme with the Monte-Carlo scheme, the per- 
centage errors for the two schemes were computed. Results are shown in Figure 5. 
It can be observed that both schemes overpredict fluid temperature while they un- 
derpredict solid temperature, and the difference between the two schemes is small. 
However, from the formulation, it is clear that this scheme requires fewer number 
of computations than the Monte-Carlo scheme. 

CONCLUSION 

A deterministic and a probabilistic one-dimensional finite element model for heat 
transfer processes in porous heat exchangers has been presented. A set of numeri- 
cal experiments have been performed to validate the model. This formulation leads 
to an asymmetric global coefficient matrix. Numerical experiments show that this 
scheme agrees well with the analytical and the Monte-Carlo methods. However, for 
mean and standard deviations, this scheme requires fewer number of computations 
in comparison to the Monte- Carlo scheme. 
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Figure 1. One dimensional porous heat exchanger 
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Figure 2. Temperature distribution 
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Figure 3. Effect of k s on temperature distribution t 
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Figure 5. Error in mean values for temperatures 
of Monte-Carlo and this scheme 
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